
library(ape)
library(plotrix)


i <- 1

phylogeny <- read.table(paste("sweep",i,"/sweep",i,".phylogeny.log",sep=""),header=FALSE)

tree <- read.tree(text=as.character(phylogeny[nrow(phylogeny),3]))


# Read the patterns
patterns <- read.csv(paste("sweep",i,"/sweep",i,".patterns",sep=""),header=FALSE)

tips <- as.numeric(tree$tip.label)

phylopatterns <- patterns[patterns[,1]%in%tips,]


# PCA them
pat <- phylopatterns[,2:245]
# PCA
pat_pca <- prcomp(pat)
red_col <- predict(pat_pca)[,1]
grn_col <- predict(pat_pca)[,2]
blu_col <- predict(pat_pca)[,3]
red_col <- red_col+abs(min(red_col))
red_col <- red_col/max(red_col)
grn_col <- grn_col+abs(min(grn_col))
grn_col <- grn_col/max(grn_col)
blu_col <- blu_col+abs(min(blu_col))
blu_col <- blu_col/max(blu_col)

# Assign colors to edges
edge.colors <- NULL
j <- 1
for(i in 1:nrow(tree$edge)){
        if(tree$edge[i,2]==j){
		ind <- which(phylopatterns[,1]==as.numeric(tree$tip.label[j]))
                edge.colors[i] <- rgb(red_col[ind],grn_col[ind],blu_col[ind])
                j <- j+1
        } else {
		ind <- which(phylopatterns[,1]==as.numeric(tree$tip.label[j]))
                edge.colors[i] <- rgb(red_col[ind],grn_col[ind],blu_col[ind])
        }
}

i <- 1
pdf(paste("sweep",i,"/sweep",i,".phylogeny.pdf",sep=""),width=2,height=2,pointsize=8)
par(mar=c(3,0,0,0.2))
plot(tree,type="c",x.lim=c(0,7300),no.margin=FALSE,root.edge=TRUE,show.node.label=TRUE,edge.color=edge.colors,show.tip.label=FALSE)
par(mgp=c(2,1,0))
axis(side=1,labels=c(0,5,10,15,20),at=c(0,1820,3640,5460,7300))
mtext(side=1,line=2,"Years")
dev.off()


i <- 1
        crypts <- read.csv(paste("sweep",i,"/sweep",i,".crypts",sep=""))
	# A clone is a unique ms pattern
        clones <- read.csv(paste("sweep",i,"/sweep",i,".grid",sep=""))
	# pick a random color for each clone
	#clonecolors <- NULL
	#clonecount <- max(clones[,2])
	ccc <- NULL
	for(j in 1:nrow(clones)){
	   if(clones[j,3]!=-1){
		ind <- which(clones[j,3]==phylopatterns[,1])
	   	ccc[j] <- rgb(red_col[ind],grn_col[ind],blu_col[ind]) 
	   } else {
		ccc[j] <- rgb(0,0,0) 
	   }
	}
	clones <- cbind(clones,ccc)
	#clonecolors <- cbind(c(1:clonecount),red_col,grn_col,blu_col)
	# plot the matrix

	mat <- matrix(crypts[,1],nrow=256,ncol=256,byrow=TRUE)
	matcol <- mat
	for(x in 1:256){
	for(y in 1:256){
	if(x==256&y==256){
		matcol[x,y]<-rgb(255,255,255,max=255)
	} else {
matcol[x,y]<-as.character(clones[(x-1)*256+y,6])
	}
	}
	}

i <- 1
pdf(paste("sweep",i,"/sweep",i,".grid.pdf",sep=""),width=3,height=3,pointsize=8)
color2D.matplot(mat,cellcolors=matcol,xlab=NA,ylab=NA,do.hex=TRUE,border=NA,axes=FALSE)
dev.off()




#matcol[x,y]<-rgb(clonecolors[clones[(x-1)*256+y,2],2],clonecolors[clones[(x-1)*256+y,2],3],clonecolors[clones[(x-1)*256+y,2],4])
